function res = h2constrel(u, udt, udtdt, params)
    A = params.A;
    Nx = params.Nx;
    Ny = params.Ny;
    Nt = Nx.*Ny;

    T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
    P1 = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);
    P2 = reshape(u(4*Nt+1 : 5*Nt), Nx, Ny);


%     [PP1, PP2] = h2icgen(T, ux, uy, params);
% 
%     DP1 = (PP1);
%     DP2 = (PP2);
%     ep1=sqrt(sum(sum(DP1.^2)));
%     ep2=sqrt(sum(sum(DP2.^2)));
%     disp(sprintf('%f %f %f %f', ep1, ep2, sqrt(sum(sum(P1.^2))), sqrt(sum(sum(P2.^2))) ));
% 
% 
%     P1mean = mean(reshape(P1./PP1, 1, params.Nx.*params.Ny));
%     P1dev = std(reshape(P1./PP1, 1, params.Nx.*params.Ny));
% 
%     P2mean = mean(reshape(P2./PP2, 1, params.Nx.*params.Ny));
%     P2dev = std(reshape(P2./PP2, 1, params.Nx.*params.Ny));
% 
%     res = [P1mean P1dev P2mean P2dev];
% 
% return;
    Tdt  = reshape(udt(0*Nt+1 : 1*Nt), Nx, Ny);
    uxdt = reshape(udt(1*Nt+1 : 2*Nt), Nx, Ny);
    uydt = reshape(udt(2*Nt+1 : 3*Nt), Nx, Ny);
    P1dt = reshape(udt(3*Nt+1 : 4*Nt), Nx, Ny);
    P2dt = reshape(udt(4*Nt+1 : 5*Nt), Nx, Ny);

    Tdtdt  = reshape(udtdt(0*Nt+1 : 1*Nt), Nx, Ny);
    uxdtdt = reshape(udtdt(1*Nt+1 : 2*Nt), Nx, Ny);
    uydtdt = reshape(udtdt(2*Nt+1 : 3*Nt), Nx, Ny);
    P1dtdt = reshape(udtdt(3*Nt+1 : 4*Nt), Nx, Ny);
    P2dtdt = reshape(udtdt(4*Nt+1 : 5*Nt), Nx, Ny);

    ut = sqrt(1+ux.^2+uy.^2);
    utdt = (uxdt.*ux + uydt.*uy)./ut;
    utdtdt = (1+ux.^2+uy.^2).^(-3/2).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+(1+ux.^2+uy.^2).*(ux.*uxdtdt+uy.*uydtdt));


    Tdx  = dx(T);
    Tdy  = dy(T);

    utdx = dx(ut);
    utdy = dy(ut);

    uxdx = dx(ux);
    uxdy = dy(ux);

    uydx = dx(uy);
    uydy = dy(uy);

    % P1dx = dx(P1);
    % P1dy = dy(P1);

    % P2dx = dx(P2);
    % P2dy = dy(P2);


    Tdtdx  = dx(Tdt);
    Tdtdy  = dy(Tdt);
    Tdxdx  = dx(Tdx);
    Tdxdy  = dy(Tdx);
    Tdydy  = dy(Tdy);

    utdtdx  = dx(utdt);
    utdtdy  = dy(utdt);
    utdxdx  = dx(utdx);
    utdxdy  = dy(utdx);
    utdydy  = dy(utdy);

    uxdtdx  = dx(uxdt);
    uxdtdy  = dy(uxdt);
    uxdxdx  = dx(uxdx);
    uxdxdy  = dy(uxdx);
    uxdydy  = dy(uxdy);

    uydtdx  = dx(uydt);
    uydtdy  = dy(uydt);
    uydxdx  = dx(uydx);
    uydxdy  = dy(uydx);
    uydydy  = dy(uydy);

    % P1dtdx  = dx(P1dt);
    % P1dtdy  = dy(P1dt);
    % P1dxdx  = dx(P1dx);
    % P1dxdy  = dy(P1dx);
    % P1dydy  = dy(P1dy);

    % P2dtdx  = dx(P2dt);
    % P2dtdy  = dy(P2dt);
    % P2dxdx  = dx(P2dx);
    % P2dxdy  = dy(P2dx);
    % P2dydy  = dy(P2dy);


    PP1 = (1/8).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(4.*uxdtdx+4.*A.*uxdtdx+(-4).*(2+A).*ux.^4.*uxdt.*uxdy.*uy.^3+4.*(2+A).*uxdt.*uxdy.*uy.^7+2.*(2+A).*(uxdt.^2+uxdy.^2).*uy.^6.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*ux.^3.*uy.*(7.*uxdtdy+(-9).*uydtdx)+2.*(1+A).*ux.*uy.^3.*(9.*uxdtdy+(-7).*uydtdx)+2.*(1+A).*ux.^5.*uy.*(3.*uxdtdy+(-5).*uydtdx)+2.*(1+A).*ux.*uy.^5.*(5.*uxdtdy+(-3).*uydtdx)+8.*(1+A).*ux.*uy.*(uxdtdy+(-1).*uydtdx)+16.*(1+A).*ux.^3.*uy.^3.*(uxdtdy+(-1).*uydtdx)+(-4).*uydtdy+(-4).*A.*uydtdy+(-2).*(1+A).*ux.^6.*((-2).*uxdtdx+uydtdy)+(-26).*(1+A).*ux.^2.*uy.^2.*((-1).*uxdtdx+uydtdy)+(-2).*(1+A).*uy.^6.*((-1).*uxdtdx+2.*uydtdy)+(-2).*(1+A).*ux.^4.*((-7).*uxdtdx+4.*uydtdy)+(-2).*(1+A).*ux.^2.*((-7).*uxdtdx+5.*uydtdy)+(-2).*(1+A).*ux.^4.*uy.^2.*((-7).*uxdtdx+6.*uydtdy)+(-2).*(1+A).*ux.^2.*uy.^4.*((-6).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^2.*((-5).*uxdtdx+7.*uydtdy)+(-2).*(1+A).*uy.^4.*((-4).*uxdtdx+7.*uydtdy)+(-4).*(2+A).*ux.^7.*uydt.*uydx+4.*(2+A).*ux.^3.*uy.^4.*uydt.*uydx+(-4).*ux.*(T.*uxdt+(-2).*(1+A).*uxdt.*uxdx+(3+A).*uydt.*uydx)+(-2).*ux.^5.*(T.*uxdt+(-1).*(1+A).*uxdt.*uxdx+(15+7.*A).*uydt.*uydx)+ux.^3.*((-6).*T.*uxdt+8.*(1+A).*uxdt.*uxdx+(-2).*(17+7.*A).*uydt.*uydx)+(-2).*(2+A).*ux.^6.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdtdt+3.*uxdxdx+2.*uxdydy+(-5).*uydxdy)+2.*(1+A).*ux.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(3.*uxdtdt+uxdxdx+2.*uxdydy+(-3).*uydxdy)+4.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+uxdxdx+uxdydy+(-2).*uydxdy)+(-4).*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-6).*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+(-2).*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdtdt+(-1).*uxdxdx+uydxdy)+4.*(2+A).*ux.^6.*uy.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+(-4).*(2+A).*ux.^2.*uy.^5.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdy+uydx).*((-1).*uxdx+uydy)+4.*(2+A).*ux.^5.*uy.^2.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-4).*(2+A).*ux.*uy.^6.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+(-2).*ux.*uy.^2.*(5.*T.*uxdt+(-13).*uxdt.*uxdx+(-11).*A.*uxdt.*uxdx+uxdy.*uydt+(-1).*A.*uxdy.*uydt+6.*uydt.*uydx+2.*A.*uydt.*uydx+(-1).*((-1)+A).*uxdt.*uydy)+ux.^3.*uy.^2.*((-8).*T.*uxdt+8.*uxdt.*uxdx+8.*A.*uxdt.*uxdx+8.*uxdy.*uydt+8.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+8.*(1+A).*uxdt.*uydy)+(-2).*ux.*uy.^4.*(3.*T.*uxdt+(-13).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+6.*uxdy.*uydt+2.*A.*uxdy.*uydt+uydt.*uydx+A.*uydt.*uydx+2.*(3+A).*uxdt.*uydy)+(-4).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+(-2).*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdx.^2+uxdy.^2+(-1).*uydt.^2+2.*uxdy.*uydx+2.*uxdx.*uydy+(-1).*uydy.^2)+(-2).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+(-2).*uxdy.*uydx+(-1).*uydx.^2+(-2).*uxdx.*uydy+uydy.^2)+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+8.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+4.*(3+A).*uxdy.^2+(-1).*uydt.^2+(-1).*A.*uydt.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdt.^2+2.*A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-2).*A.*uxdy.^2+(-2).*uydt.^2+(-2).*A.*uydt.^2+2.*A.*uydx.^2+(-1).*uydy.^2+(-1).*A.*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+A.*uxdt.^2+uxdx.^2+A.*uxdx.^2+(-14).*uydt.^2+(-8).*A.*uydt.^2+(-12).*uydx.^2+(-4).*A.*uydx.^2+(-1).*(1+A).*uydy.^2+2.*T.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+6.*uxdx.^2+6.*A.*uxdx.^2+(-5).*uydt.^2+(-3).*A.*uydt.^2+(-6).*(1+A).*uydy.^2+4.*T.*((-1).*uxdx+uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(14.*uxdt.^2+10.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*((-4)+A).*uxdy.^2+(-4).*uydt.^2+(-4).*A.*uydt.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+4.*A.*uxdt.^2+3.*uxdx.^2+3.*A.*uxdx.^2+(-2).*A.*uxdy.^2+(-14).*uydt.^2+(-10).*A.*uydt.^2+(-8).*uydx.^2+2.*A.*uydx.^2+(-3).*uydy.^2+(-3).*A.*uydy.^2+6.*T.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+2.*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+(-2).*T.*(uxdy+(-1).*uydx)+3.*uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-9).*uydx.*uydy+(-5).*A.*uydx.*uydy+(1+A).*uxdy.*(uxdx+uydy))+(-2).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.*uydt+3.*A.*uxdt.*uydt+2.*T.*(uxdy+(-1).*uydx)+uxdx.*uydx+A.*uxdx.*uydx+uydx.*uydy+A.*uydx.*uydy+(-1).*uxdy.*((9+5.*A).*uxdx+((-3)+A).*uydy))+uy.^3.*(2.*(17+7.*A).*uxdt.*uxdy+2.*uydt.*(3.*T+(-4).*(1+A).*uydy))+4.*uy.*((3+A).*uxdt.*uxdy+uydt.*(T+(-2).*(1+A).*uydy))+2.*uy.^5.*((15+7.*A).*uxdt.*uxdy+uydt.*(T+(-1).*(1+A).*uydy))+(-2).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(2.*T.*(uxdy+(-1).*uydx)+(-1).*uxdy.*((5+3.*A).*uxdx+(1+3.*A).*uydy)+uydx.*((1+3.*A).*uxdx+(5+3.*A).*uydy))+2.*ux.^4.*uy.*((1+A).*uxdt.*uxdy+2.*(3+A).*uxdt.*uydx+uydt.*(3.*T+2.*(3+A).*uxdx+(-1).*(13+9.*A).*uydy))+2.*ux.^2.*uy.*(2.*(3+A).*uxdt.*uxdy+(-1).*((-1)+A).*uxdt.*uydx+uydt.*(5.*T+(-1).*((-1)+A).*uxdx+(-1).*(13+11.*A).*uydy))+4.*ux.^2.*uy.^3.*((3+A).*uxdt.*uxdy+(-2).*((1+A).*uxdt.*uydx+uydt.*((-1).*T+(1+A).*(uxdx+uydy))))+(-4).*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-6).*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-2).*(1+A).*uy.^5.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdxdy+uydtdt+uydydy)+(-4).*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*((-2).*uxdxdy+2.*uydtdt+uydxdx+uydydy)+(-2).*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-3).*uxdxdy+3.*uydtdt+2.*uydxdx+uydydy)+(-2).*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*((-5).*uxdxdy+5.*uydtdt+2.*uydxdx+3.*uydydy));

    PP2 = (1/4).*T.*(1+ux.^2+uy.^2).^(-3/2).*(2+ux.^2+uy.^2).^(-1).*(2.*uxdtdy+2.*A.*uxdtdy+2.*(1+A).*ux.*uxdtdx.*uy.^5+4.*(1+A).*uxdtdy.*uy.^6+(-4).*(2+A).*ux.*uxdt.*uxdy.*uy.^6+2.*(1+A).*ux.*uxdxdy.*uy.^4.*(1+ux.^2+uy.^2).^(1/2)+(-2).*(2+A).*ux.*(uxdt.^2+uxdy.^2).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*(1+A).*(uxdtdt+uxdydy).*uy.^5.*(1+ux.^2+uy.^2).^(1/2)+2.*uydtdx+2.*A.*uydtdx+4.*(1+A).*ux.^6.*uydtdx+10.*(1+A).*ux.^2.*uy.^2.*(uxdtdy+uydtdx)+4.*(1+A).*uy.^2.*(2.*uxdtdy+uydtdx)+2.*(1+A).*ux.^2.*uy.^4.*(3.*uxdtdy+uydtdx)+2.*(1+A).*uy.^4.*(5.*uxdtdy+uydtdx)+4.*(1+A).*ux.^2.*(uxdtdy+2.*uydtdx)+2.*(1+A).*ux.^4.*uy.^2.*(uxdtdy+3.*uydtdx)+2.*(1+A).*ux.^4.*(uxdtdy+5.*uydtdx)+2.*(1+A).*ux.^5.*uy.*uydtdy+2.*(1+A).*ux.*uy.*(uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.^3.*(uxdtdx+uydtdy)+2.*(1+A).*ux.*uy.^3.*(2.*uxdtdx+uydtdy)+2.*(1+A).*ux.^3.*uy.*(uxdtdx+2.*uydtdy)+(-4).*(2+A).*ux.^6.*uy.*uydt.*uydx+(-2).*(2+A).*ux.^5.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uydt.^2+uydx.^2)+2.*(1+A).*ux.^5.*(1+ux.^2+uy.^2).^(1/2).*(uydtdt+uydxdx)+2.*(1+A).*ux.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdxdy+uydtdt+uydxdx)+2.*(1+A).*ux.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdxdy+2.*(uydtdt+uydxdx))+2.*(1+A).*ux.^4.*uy.*(1+ux.^2+uy.^2).^(1/2).*uydxdy+2.*(1+A).*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+uydxdy)+2.*(1+A).*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(2.*uxdtdt+2.*uxdydy+uydxdy)+2.*(1+A).*ux.^2.*uy.*(1+ux.^2+uy.^2).^(1/2).*(uxdtdt+uxdydy+2.*uydxdy)+4.*(2+A).*ux.^5.*uy.^2.*(uxdt.*uydx+uydt.*(uxdx+(-1).*uydy))+4.*(2+A).*ux.^2.*uy.^5.*((-1).*uxdt.*uxdx+uxdy.*uydt+uxdt.*uydy)+4.*(2+A).*ux.^4.*uy.^3.*((-1).*uxdt.*uxdx+uydt.*(uxdy+(-1).*uydx)+uxdt.*uydy)+ux.^4.*uy.*((-4).*(2+A).*uxdt.*uxdx+uydt.*(2.*(3+A).*uxdy+(-1).*(19+9.*A).*uydx)+2.*(3+A).*uxdt.*uydy)+uy.^5.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(5+3.*A).*uxdt.*uydy)+uy.*((-2).*T.*uxdt+(-1).*uxdt.*uxdx+A.*uxdt.*uxdx+5.*uxdy.*uydt+3.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+(5+3.*A).*uxdt.*uydy)+uy.^3.*((-4).*T.*uxdt+(-2).*uxdt.*uxdx+2.*A.*uxdt.*uxdx+10.*uxdy.*uydt+6.*A.*uxdy.*uydt+(-1).*uydt.*uydx+A.*uydt.*uydx+2.*(5+3.*A).*uxdt.*uydy)+ux.^2.*uy.*((-2).*T.*uxdt+(-11).*uxdt.*uxdx+(-5).*A.*uxdt.*uxdx+11.*uxdy.*uydt+5.*A.*uxdy.*uydt+(-12).*uydt.*uydx+(-4).*A.*uydt.*uydx+(11+5.*A).*uxdt.*uydy)+ux.^2.*uy.^3.*((-2).*T.*uxdt+(-19).*uxdt.*uxdx+(-9).*A.*uxdt.*uxdx+21.*uxdy.*uydt+11.*A.*uxdy.*uydt+(-9).*uydt.*uydx+(-3).*A.*uydt.*uydx+(21+11.*A).*uxdt.*uydy)+(-4).*(2+A).*ux.^4.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*((-1).*uxdt.*uydt+(-1).*uxdx.*uydx+uydx.*uydy)+ux.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uydx+(5+3.*A).*(uxdt.*uydt+uxdx.*uydx)+((-1)+A).*uydx.*uydy)+(-2).*(2+A).*ux.^3.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.^2+uxdx.^2+uydt.^2+(-2).*uxdy.*uydx+(-2).*uxdx.*uydy+uydy.^2)+(-1).*ux.*uy.*(1+ux.^2+uy.^2).^(1/2).*(5.*uxdt.^2+3.*A.*uxdt.^2+2.*uxdx.^2+2.*uxdy.^2+5.*uydt.^2+3.*A.*uydt.^2+(-4).*uxdy.*uydx+2.*uydx.^2+(-4).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.^3.*uy.*(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.^2+2.*A.*uxdt.^2+2.*(2+A).*uxdx.^2+9.*uydt.^2+5.*A.*uydt.^2+(-6).*uxdy.*uydx+(-2).*A.*uxdy.*uydx+6.*uydx.^2+2.*A.*uydx.^2+(-2).*(3+A).*uxdx.*uydy+2.*uydy.^2)+(-1).*ux.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.^2+5.*A.*uxdt.^2+2.*uxdx.^2+2.*(3+A).*uxdy.^2+4.*uydt.^2+2.*A.*uydt.^2+(-2).*(3+A).*uxdy.*uydx+(-2).*(3+A).*uxdx.*uydy+2.*(2+A).*uydy.^2)+4.*(2+A).*ux.^2.*uy.^4.*(1+ux.^2+uy.^2).^(1/2).*(uxdt.*uydt+uxdy.*((-1).*uxdx+uydy))+ux.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+6.*uxdx.*uydx+2.*A.*uxdx.*uydx+(-2).*T.*(uxdy+2.*uydx)+4.*A.*uydx.*uydy+(-1).*((-1)+A).*uxdy.*((-1).*uxdx+uydy))+ux.^2.*uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(20.*uxdt.*uydt+12.*A.*uxdt.*uydt+9.*uxdx.*uydx+3.*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+(-9).*uydx.*uydy+(-3).*A.*uydx.*uydy+3.*(3+A).*uxdy.*((-1).*uxdx+uydy))+(-4).*(2+A).*ux.^3.*uy.^4.*(uxdt.*uxdy+(-1).*uxdt.*uydx+uydt.*((-1).*uxdx+uydy))+(1+ux.^2+uy.^2).^(1/2).*(4.*uxdt.*uydt+4.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(uxdy+uydx)+uydx.*uydy+3.*A.*uydx.*uydy+uxdy.*((1+3.*A).*uxdx+(-1).*((-1)+A).*uydy))+ux.^5.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+ux.*(((-1)+A).*uxdt.*uxdy+(5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy))+uy.^2.*(1+ux.^2+uy.^2).^(1/2).*(9.*uxdt.*uydt+7.*A.*uxdt.*uydt+uxdx.*uydx+(-1).*A.*uxdx.*uydx+(-2).*T.*(2.*uxdy+uydx)+(-1).*uydx.*uydy+A.*uydx.*uydy+2.*uxdy.*(2.*A.*uxdx+(3+A).*uydy))+uy.^4.*(1+ux.^2+uy.^2).^(1/2).*((-2).*T.*uxdy+(5+3.*A).*uxdt.*uydt+uxdy.*(((-1)+A).*uxdx+(5+3.*A).*uydy))+ux.^3.*uy.^2.*((-3).*(3+A).*uxdt.*uxdy+(21+11.*A).*uxdt.*uydx+uydt.*((-2).*T+(21+11.*A).*uxdx+(-1).*(19+9.*A).*uydy))+ux.*uy.^2.*((-4).*(3+A).*uxdt.*uxdy+(11+5.*A).*uxdt.*uydx+uydt.*((-2).*T+(-1).*(11+5.*A).*((-1).*uxdx+uydy)))+ux.^3.*(((-1)+A).*uxdt.*uxdy+2.*((5+3.*A).*uxdt.*uydx+uydt.*((-2).*T+(5+3.*A).*uxdx+((-1)+A).*uydy)))+ux.*uy.^4.*((-1).*(19+9.*A).*uxdt.*uxdy+2.*((3+A).*uxdt.*uydx+uydt.*((3+A).*uxdx+(-2).*(2+A).*uydy))));

    DP1 = (PP1);
    DP2 = (PP2);
    ep1=sqrt(sum(sum(DP1.^2)));
    ep2=sqrt(sum(sum(DP2.^2)));
    disp(sprintf('%f %f %f %f', ep1, ep2, sqrt(sum(sum(P1.^2))), sqrt(sum(sum(P2.^2))) ));


    P1mean = mean(reshape(P1./PP1, 1, params.Nx.*params.Ny));
    P1dev = std(reshape(P1./PP1, 1, params.Nx.*params.Ny));

    P2mean = mean(reshape(P2./PP2, 1, params.Nx.*params.Ny));
    P2dev = std(reshape(P2./PP2, 1, params.Nx.*params.Ny));

    res = [P1mean P1dev P2mean P2dev];
end
